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-Abstract- 

We present two quantitative behavioral equivalences over species of a chemical reaction network 
(CRN) with semantics based on ordinary differential equations. Forward CRN bisimulation iden¬ 
tifies a partition where each equivalence class represents the exact sum of the concentrations 
of the species belonging to that class. Backward CRN bisimulation relates species that have 
identical solutions at all time points when starting from the same initial conditions. Both no¬ 
tions can be checked using only CRN syntactical information, i.e., by inspection of the set of 
reactions. We provide a unified algorithm that computes the coarsest refinement up to our bisim¬ 
ulations in polynomial time. Further, we give algorithms to compute quotient CRNs induced 
by a bisimulation. As an application, we find significant reductions in a number of models of 
biological processes from the literature. In two cases we allow the analysis of benchmark models 
which would be otherwise intractable due to their memory requirements. 

1998 ACM Subject Classification D.2.4, G.1.7, J.2 

Keywords and phrases Chemical reaction networks - ordinary differential equations - bisimula¬ 
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1 Introduction 

At the interface between computer science and systems biology is the idea that biological 
systems can be interpreted as computational processes I22HH], leading to a number of formal 
methods applied to study biomolecular systems HEW In this context, chemical reaction 
networks (CRNs), a popular mathematical model of interaction in natural sciences, can also 
be seen as a kernel concurrent language for natural programming. 

In this paper we present, for the first time to our knowledge, quantitative bisimulation 
equivalences for CRNs with the well-known interpretation based on ordinary differential 
equations (ODEs). (To make the paper self-contained, all background is given in Section !) 
In this semantics, each species is associated with an ODE giving the deterministic evolution 
of its concentration starting from an initial condition. Our bisimulations are equivalences 
over species that induce a reduced CRN that exactly preserves the dynamics of the original 
one. This is an important goal, especially in order to cope with the potentially very large 
number of species and reactions in many biological networks cans]. 

We study two equivalences, developed in the Larsen-Skou style of probabilistic bisim¬ 
ulation EH], that are based on two distinct ideas of observable behavior. Forward CRN 
bisimulation yields an aggregated ODE where the solution gives the exact sum of the concen¬ 
trations of the species belonging to each equivalence class. In backward CRN bisimulation, 
instead, equivalent species have the same solution at all time points; in other words, backward 
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Forward and Backward Bisimulations for Chemical Reaction Networks 


CRN bisimulation relates species whose ODE solutions are equal whenever they start from 
identical initial conditions. The use of “forward” and “backward” has a long tradition in 
models of computation based on labelled transition systems jl8l . In the case of quantitative 
variants, for instance those defined for process algebra with a continuous-time Markov chain 
(CTMC) semantics {251 SB! 181 IT] . forward bisimulations are equivalences that induce a CTMC 
aggregation in the sense of ordinary lumpability [7], where the probability of an equivalence 
class is equal to the sum of the probabilities of the states belonging to that class. This is found 
by checking conditions on the outgoing transitions of related states in the transition diagram. 
A backward bisimulation induces a CTMC aggregation in the sense of weak lumpability WL 
where all states in the same equivalence class have a time-invariant conditional probability 
distribution; exact lumpability is a special case where the conditional probability distribution 
is uniform, in the sense that any two states of each equivalence class have the same probability 
at any point in time whenever they have the same initial probabilities. It is found by relating 
states according to conditions on their predecessor states [HU GS| Ej . 


Despite being similar in spirit, technically our bisimulations are fundamentally different 
for two reasons. First, they concern a continuous-state semantics based on ODEs instead of 
a discrete-state CTMC. Second, they operate at the structural, syntactical level, because 
they are defined with quantities that can be computed by only inspecting the reactions 
of a CRN. Still, the repercussions of our bisimulations on the semantics are explained 
according to certain theories of aggregation. In particular, forward CRN bisimulation yields 
an aggregated system in the sense of ODE lumpability [mug. This theory covers linear 
transformations of the original state variables in general; here we consider an instance, which 
we call ordinary fluid lumpability , where the transformation is induced by a partition of 
state variables. (Forward bisimulation is presented in Section 3.1 ) Backward bisimulation 
(presented in Section 3.2) is related to exact fluid lumpability , introduced in the context 
of process algebra with fluid semantics [3B], identifying process terms with the same ODE 
solution when initialized equally. The disadvantage of forward CRN bisimulation is that it is 
lossy (yet exact) because, similarly to the forward stochastic analogues, from the aggregated 
ODE system in general it is not possible to recover the solutions for the individual species 
within the same equivalence class. On the other hand, it does not impose restrictions on the 
initial conditions, which instead are present in our backward variant. As a further important 
difference, forward CRN bisinrulation (again, like its stochastic analogues) turns out to be 
a sufficient condition for ODE lumpability. Instead, backward CRN bisimulation enjoys a 
full characterization, in the sense that there exists a backward CRN bisinrulation between 
two species if and only if they have the same ODE solutions (provided that they start from 
equal initial conditions). More in general, by means of a number of examples we will show 
that the two equivalences are complementary because not comparable. In other words, there 
exist models that can be reduced up to forward CRN bisimulation but not by the backward 
variant, and vice versa; at the same time, there are models that can be reduced by both. 


To enhance the usefulness of these notions, we present (in Section [5| a template partition- 
refinement algorithm that is parametric with respect to the equivalence of interest, computing 
the coarsest refinement up to either variant in polynomial time. To use our equivalences as 
an automatic model reduction tool, we further give two algorithms (in Section [4]) that provide 
the quotient CRN induced by either bisimulation. With a prototype implementation available 
at http://sysma.imtlucca.it/crnreducer/, we show (in Section [6]) that we are able to 
reduce a number of case studies taken from the literature. Our bisimulations yielded quotient 
CRNs with number of reactions and species up to four orders of magnitude smaller than the 
original CRNs, leading to speed-ups in the ODE solution runtimes of up to five orders of 
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magnitude. In two cases, it was possible to analyze models that were otherwise intractable 
directly within our experimental environment due to excessive memory requirements. 

Related work. Behavioral equivalences have been recently proposed in m for comparing 
CRNs; however, the analysis is carried out at the qualitative level, i.e., ignoring the dynamical 
evolution. In [36] is introduced the notion of label equivalence for process algebra with fluid 
semantics, which captures exact fluid lumpability (processes are equivalent whenever their 
ODE solutions are equal at all time points). However, unlike backward CRN bisimulation, 
label equivalence is only a sufficient condition for ODE reduction. Indeed, it works at a 
coarser level of granularity as it relates sets of ODE variables, each corresponding to the 
behavior of a sequential process. Instead, backward CRN bisimulation relates individual 
ODE variables. Further, the conditions for equivalence, specific to the process algebra, are 
difficult to check automatically because of the universal quantifiers over the ODE variables. 
More important, no algorithm for computing the coarsest partition was developed. Similar 
considerations apply to the process-algebra specific ordinary fluid lumpability in m- 

Cardelli’s notion of emulation between two CRNs is a (structural) mapping of species and 
reactions that, like backward CRN bisimulation, guarantees the equality between the ODE 
solutions at all time points m- An emulation requires a source and a target CRN — the 
modeler is intended to have the suspicion that, for some given CRN, another CRN might be 
related to it. But emulation cannot be used when one wants to discover equivalences between 
species within the same given CRN. Thus, emulation is not useful for model reduction because 
a-priori information about the structure of a quotient CRN is not available. Furthermore, 
no algorithm is provided in m to find emulations automatically. Since backward CRN 
bisimulation fully characterizes exact fluid lumpability, it is not difficult to show that backward 
CRN bisimulation generalizes emulation in the sense that any emulation between two CRNs 
can be understood in terms of a backward CRN bisinrulation over the species of a “union 
CRN” that contains all the reactions of the two CRNs of interest. 

Model reductions have been studied in related models for biomolecular networks (e.g. urn 
HD US), most notably for rule-based systems such as BioNetGen [5] and the n calculus BZI- 
These offer an intensional modeling approach, by providing graph-rewrite rules of interaction 
instead of a complete enumeration of all chemical reactions involved. Differential fragments 
for k are self-consistent aggregates found by a static analysis on the model, identifying sums 
of chemical species for which an ODE system can be explicitly written m- In this sense, they 
are analogous to our CRN bisimulations, but with notable differences. First, fragmentation 
works directly at the rule-based level. This has the advantage that the analysis is performed on 
a set of rewrite rules, which is typically much more compact than the fully enumerated CRN. 
However, fragmentation is domain-specific, hence the model must be conveniently expressed 
as a biomolecular system (e.g., with complex formation or internal state modification). On 
the other hand, CRN bisimulations work for a generic language-independent CRN, which 
however must be explicitly given. Further, unlike CRN bisimulations, fragmentation is 
performed on a “static” view of the model, without information on the reaction rates. The 
ODE aggregations of both forward CRN bisimulation and fragmentation introduce loss of 
information (in contrast to backward CRN bisimulation). But, unlike our forward variant, in 
fragmentation the same species may be present in more than one fragment. Additionally, 
species may occur in fragments with multiplicity numbers. Thus, fragmentation can be seen 
as a form of improper lumping that is not necessarily induced by a partition of the original 
state-space variables m ■ Overall, because of these differences, it is not difficult to find 
models that can be reduced by our CRN bisimulations and not by fragmentation, and vice 
versa. This is presented in detail in Section [6| 
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2 Background 

Notation. We write A —>■ B and B A for the functions from A to B. When f £ A —>• B and 
a £ A, we set f a := /(a). Moreover, for any X C A and b £ B, we define f(X) := {b £ B \ 
3 a £ X.(f(a ) = b)}. Sets and multisets are denoted by {...} and {| ... |}, respectively. Also, 
we shall not distinguish among an equivalence relation and the partition induced by it, and 
shall use the symbol to denote the equivalence relation with H = Finally, given 

two partitions Hi and H 2 of a given set S , we say that H 1 is a refinement of H 2 if for any 
Hi £ Hi there exists a (unique) H 2 £ H 2 such that Hi C H 2 . 

2.1 Chemical Reaction Networks 

Formally, a CRN (S, R) is a pair consisting of a finite set of species S taken from a countable 
infinite universe of all species, and a finite set of chemical reactions R. A reaction is a triple 
written in the form p -4- n, where p and 7 r are the multisets of species reactants and products , 
respectively, and a > 0 is the reaction rate. In particular, we focus on basic chemistry where 
only elementary reactions are considered, where at most two reactants (possibly of the same 
species) interact. No restrictions are instead imposed on products. Several models found 
in the literature (including those discussed in Section [ 6 | belong to this class. Also, this is 
consistent with the physical considerations which stipulate that reactions with more than 
two reactants are very unlikely to occur in nature m■ We denote by p(X) the multiplicity 
of species X in the multiset p , and by A4S(S) the set of finite multisets of species in S. To 
adhere to standard chemical notation, we shall use the operator + to denote multiset union, 
e.g., X + Y + Y (or just X + 2Y) denotes the multiset {|AT, Y, F|}. We may also use X to 
denote either the species X or the singleton {] X |}. 

The (autonomous) ODE system V = F(V) underlying a CRN (S,R) is F : R > 0 —► R s , 
where each component Fx, with X £ S is defined as: 

p x(V) := £ MX) - p(X)) a JJ V y p(Y) . 

p-^AireR Yes 

This represents the well-known mass-action kinetics, where the reaction rate is proportional 
to the concentrations of the reactants involved. Since the ODE system of a CRN is given by 
polynomials, the vector field F is locally Lipschitz. Hence, the theorem of Picard-Lindelof 
ensures that for any H( 0 ) £ R > 0 there exists a unique non-continuable solution of V = F(V). 

► Example 1. We now provide a simple CRN, ( S e ,R e ), with S e = {A, B,C, D, E} and 
R e = {A -^4 E, B -^4 D, A+B -h> C, C+D -4 2 C+D, E+D 2E+D}, which will be 
used as a running example throughout the paper. Its ODE system is given by 

Va=-6Va-2V a Vb Vb = -6Vb-2V a Vb V c = 2V a V b + 5V c V d 
V d = 6 V b V e = 6Va + 5V e V d 

In the following, we shall assume that the universe of all species is well-ordered with 
respect to C. We then say that a function p, : S —> S is a choice function of a partition H 
of S, if p(X) = mine H for all H £ H and X £ H. Also, choice functions can be trivially 
lifted to multisets applying them element-wise, e.g., p(X + Y) = p(X) + p(Y). 

2.2 Fluid Lumpability 

Ordinary Fluid Lumpability. We start by defining the notion of ordinary fluid lumpability, 
which is an instance of ordinary lumpability for ODEs m specialized to CRNs. 
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► Definition 2 (Ordinary fluid lumpability). Let (S,R) be a CRN, F be its vector field, 
and H = {Hi ,..., H m } a partition of S. Then, H is ordinary fluid lumpable if for 
all H £ H there exists a polynomial pH in \H\ variables such that J^xeH Bx(V) = 
i Ph{^2,x^h 1 Vx, ■ ■ ■, ThxeH m Vx) for all V £ R> 0 - 

Informally, a partition H is ordinary fluid lumpable if, for each H £ H, the polynomial 
Fx(V) in the variables {Vx \ X £ S} can be rewritten into a polynomial pH in the 
variables {^2xeH I H £ R}- In particular, if H is known to be an ordinary fluid lumpable 
partition of ( S,R ) and V denotes the solution of V = F{V) subject to R(0) £ R> 0 , the 
solution of the aggregated ODE system (Wh, , ■ ■ •, W Em ) = (pH 1 {W ),..., pH m (W)) with 
HV-/(0) = J^xeH Vx(0) is such that Wnifl) = J^xgh ErW for all t £ domain (E). 

► Example 3. Consider the ODEs of (S e , R e ) of Example]!] and let Ho = {{^4}, {B}, {C, E}, 
{D}}. By applying a variable renaming consistent with the blocks of Ho, i-e., Vce = Vc+Ve, 
and by exploiting the linearity of the differential operator we get 

Va = -6V a -2V a Vb Vb = -6Vb-2V a Vb V C e = 2V a Vb+6Va+5V d Vce V d = 6V b 

That is, we obtained an ODE system in terms of block variables only. ◄ 

Exact Fluid Lumpability. We extend to CRNs the notion of exact fluid lumpability in [155] . 

► Definition 4 (Exact fluid lumpability). Let ( S,R ) be a CRN, F its vector field, and H a 
partition of S. We call V £ constant on H if Ey; = Vx^ for all H £H, and all X-, , Xj £ H. 
Then, H is exactly fluid lumpable if F(V) is constant on H whenever V is constant on H. 

► Example 5. Consider the ODEs of (S e , R e ) of Example]!] and let He = {{R, B}, {C{, {D}, 
{E}}. It is easy to see that A and B have same concentrations at all time points if initialized 
equally. In these cases, we can replace the ODEs of ( S e ,R e ) with the ones aggregated 
according to He, obtained by removing Vb and replacing all occurrences of Vb with Va- 

Va = 6 V A 2 V A V A V c = 2V a Va + 5V c V d V d = 6V a V e = 6V a + 5V e V d 

That is, we obtained a (lossless) aggregated ODE system written in terms of a variable per 
block, chosen according to C. ◄ 

We remark that the above definition expresses exact fluid lumpability in terms of properties 
of the ODE vector field of a CRN. Instead, in [55] exactly fluid lumpability was defined 
directly in terms of the desired dynamical property, i.e., that the ODE solutions within any 
equivalence class be equal at all time points. The following result is a new contribution 
showing that this dynamical property is fully characterized by the vector-field based definition. 

► Theorem 6. Let ( S,R ) be a CRN and F its vector field. A partition H of S is exactly fluid 
lumpable if and only if, for any E(0) £ M> 0 that is constant on H, the underlying solution of 
V = F(V) is such that V(t) is constant on H for all t £ domain(V). [|] 

3 CRN Bisimulations 

Both notions of fluid lumpability given in Section [2] are not convenient to be used directly 
because they involve a universal quantifier over the (uncountable) state space. We address 
this problem by providing structural conditions that concern only the reactions of a CRN. 


i 


For the sake of readability all proofs are provided in a separate appendix. 
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3.1 Forward CRN Bisimulation 


We now introduce forward CRN bisimulation, an equivalence on species that will turn out 
to induce ordinary fluid lumpability. We start with the notions of reaction and production 
rate. The former collects the rates at which the concentration of a species X decreases when 
reacting with a given partner. The latter collects the positive contribution that X exerts to 
the concentration of a species Y, again when reacting with a certain partner. 


► Definition 7 (Reaction and production rates). Let ( S,R ) be a CRN, X,Y £ S, and 
p £ A iS(S). The p-reaction rate of X , and the p-production rate of ^elements by X are 
defined respectively as 


~[X, p] (p(X) + 1) a ’ P r (^; p,Y):=(p(X) + 1) £ a ■ 7t(Y) 


x+p- 


>7 r£R 


X+p- 


> 7 t£R 


Finally, for H C S we define pr[X, p, FT] := J2yfh P r (^b P> V). 

► Definition 8 (Forward CRN Bisimulation). Let ( S,R ) be a CRN, TZ an equivalence relation 
over S and TL = S/TZ. Then, TZ is a forward CRN bisimulation (abbreviated FB) if for all 
{X, Y) £ TZ, all p £ M.S(S), and all H £ TL it holds that 


crr[X, p\ = err[y, p\ and pr[X, p, H] = pr[Y, p, H] (1) 

► Example 9. Consider TLo = {{^4}, {B}, {C, E}, {13}} of Example |3j It can be shown that 
TLo is an FB, as, e.g., crr[C,D] = err[E,D] = 5, and pr[C,D, {C,E}\ = pr[E,D, {C,E}\ = 10. 

We are interested in the coarsest FB, as well as in the coarsest one refining a given initial 
partition of species. 

► Proposition 10. Let ( S,R ) be a CRN, I a set of indices, and TZi an FB for ( S,R ), for all 
i £ I. The transitive closure of their union 7 Z= ((J ig/ TZi)* is an FB for (5, R). In particular, 
if each Ri is such that S/TZi refines some partition Q of S, then so does S/TZ. 

► Theorem 11 (Forward bisimulation implies ordinary fluid lumpability). Let ( S,R ) be a CRN. 

Then, TL is an ordinarily fluid lumpable partition of S if TL is an FB of S. 

FB is only a sufficient condition for lumpability, as discussed in the next example. (However, 
Section [6] shows that FB can be effectively applied to interesting existing models.) 

► Example 12. Consider the CRN ({F, G}, {F —4 G, G —b F}), having ODEs 

Vf = — cti Vf + Oi2 Vq Vg = — 02 Vg + Oi Vf 


If oi / 02 , TL C = {{G, F}} is not an FB, as crr[F, 0] = aq and crr[G, 0] = 02 - Nevertheless, 
the above ODE system is lumpable. Indeed, by applying the variable renaming consistent 
with TL C , i.e., Vfg = Vf + Vg, we get a single ODE for Vfg, be., Vfg = 0- ◄ 


3.2 Backward CRN Bisimulation 

We now introduce backward CRN bisimulation, an equivalence on species that will turn out 
to characterize exact fluid lumpability. We start with the notion of cumulative flux rate, 
which collects the overall contribution that reactions with a given multiset of reactants p 
exert to the concentration of a species X. 
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Figure 1 Relation among (R-reduced) CRNs and (R-lumped) semantics, with R a bisimulation. 


► Definition 13 (Cumulative flux rate). Let ( S,R ) be a CRN, X £ S, p £ A4S(S), and 
A4 C A4S(S). Then, we define 





pEM 


p 


We call fr(X, p) and fr[X, AA\ p-flux rate and cumulative A4-flux rate of X , respectively. 

► Definition 14 (Backward CRN bisimulation). Let ( S,R ) be a CRN, 7Z an equivalence 
relation over S, R = S/1Z and p the choice function of R. Then, 1Z is a backward CRN 
bisimulation (BB) if for any (A', F) £ 7Z it holds that 


fr[X, M\ = fr[F, A4] for all A4 £ {p \ p n £ R}/ 


( 2 ) 


where any two p,cr £ A4S(S) satisfy p ~-u a if p(p) = p{a). 

► Example 15. Consider Re = {{A, B}, {C}, {D}, {E}} of Example[5] We first note that 

{|A|} {|R|}, as ~ Ue relates multisets with same number of R^-equivalent species. Also, it 

can be shown that Re is a BB, as, e.g., fr[A, M] = fr[R, A4] = —6 for A4 = {{|A|}, {|A?|}}. ◄ 

As for FB, there exists a coarsest BB (that refines a given partition of S). 

► Proposition 16. Let ( S,R ) be a CRN, I a set of indices, and TZi a BB for ( S,R ), for all 
i £ I. The transitive closure of their union 1Z= (U, e / R-i)* is a BB for ( S,R ). In particular, 
if each Ri is such that S/TZi refines some partition Q of S , then so does S/TZ. 

We now state the mentioned characterization of exact fluid lumpability in terms of BB. 

► Theorem 17 (Backward bisimulation characterizes exact fluid lumpability). Let ( S,R ) be a 

CRN. Then, R is an exactly fluid lumpable partition of S if and only if R is a BB of S. 

► Remark. We wish to stress that FB and BB are not comparable: First, "Ho is not a BB, 

as fr[C,{A+I?}] = 2 and ir[E,{A+B}] =0; Second,R^is not an FB, as crr(A,R)=2 and 
err (B, B)= 0; Third, for the same reasons, {{A,B}, {C,E}, {D}} is neither an FB nor a BB. 
Similar examples on models of biological relevance are provided in Section [6] ◄ 

4 Reduced Chemical Reaction Networks up to CRN Bisimulations 

We have shown that, given a CRN and a CRN bisimulation 7Z , we can analyze the aggregated 
ODE system according to 1Z. We now provide the notion of reduced CRN from which the 
aggregated ODEs can be directly generated, as depicted in Figure |T] 

► Definition 18 (Forward reduction). Let ( S,R ) be a CRN, R an FB, and p its choice 


function. The (R, F)-reduction of (S,R) is given by (S,R)^ H ’ F ^ = R^ n ^), where 


S (H, F ) = ^( 5 ) and R (U,f) 

is defined as follows: (FI) Discard all reactions p tt such 
that p ^ p{p)\ (F2) Replace all remaining reactions p -^4 7r with p p(n); (F3) Fuse all 
reactions that have the same reactants and products by summing their rates. 
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The idea underlying forward reduction is to discard all reactions having non-representative 
reagents, and to replace the products of the remaining reactions with their representatives. 
This can be seen as a special case of Theorem 4.4 of M- 

► Example 19. Consider the FB Ho = {{A},{B},{C,E},{D}} used in Example^] The 
(Ho, -F)-reduction of ( S e ,R e ) is (with C being the representative of its block) Sx*°’ F ) = 
{A, B, C, D}, Ri Ho ’ F) = {A -^4 C, B -^4 D, A+B -4- C, C+D A 2 C+D}. Note that the 

5 

reaction E+D — >2E+D is discarded, as E is not a representative species. ◄ 

We now state that the (H, F)-reduction of an FB H induces the ODEs aggregated 
according to H. For example, the (Ho, F)-reduction of ( S e ,R e ) induces the ODEs shown in 
Example [3] if applying the renaming Vq = Vce- 

► Theorem 20 (Forward reduction induces aggregation). Let (S, R) be a CRN, H an FB and p 
its choice function. Then, (S, R) <yH ' F " 1 is computed in at most 0(\R\ ■ |5| • (log(|l?|) + log(|5|))) 
steps. Crucially, if F is the vector field of (S,R) and F the one of (S, R)^’ F \ then 

gh Bx(V) = F fJ .(Y)(J2xGH 1 Vx,--’, J2x eH m Vx) for all V £ lR>o; H £ H and Y £ H. 

For the backward reduction, the underlying idea is to keep track only of differential 
contributions that affect the representative species p(S). 

► Definition 21 (Backward reduction). Let ( S,R ) be a CRN, H a BB, and p its choice 

function. The (H, B)-reduction of (S,R) is given by (S,R)^ n,B ^ = (S^ H ’ B \R^ H ’ B " > ), where 
5 (w,b) = and R (h,b) 

is obtained as follows: (Bl) Replace all reactions p -^4 7r with 
p it where n(Xi) := -k(X f) if Xi G p(S) and 7r(Xj) := p(Xf) otherwise; (B2) Replace all 
such obtained reactions p —X 7r with p(p) p(x)\ (B3) Fuse all reactions that have the 
same reactants and products by summing their rates. 

► Example 22. Considering the CRN ( S e ,R e ) and the BB He, (Bl) changes B-^AD in 
B-^D+B, and A+B — C in A+B-^AC+B, while(B2) yields {A-^-E, A-^-D+A, A+A-^A 
C+A,C+D—^- J >-2C+D,E+D—^ J >-2E+D}. Finally, (B3) does not introduce any change. ◄ 

► Theorem 23 (Backward reduction induces aggregation). Let (S,R) be a CRN, H a BB and p 
its choice function. Then, (S,R)^’ B ^ is computed in at most 0(\R\-\S\-(\og(\R\)+\og(\S\))) 
steps. Crucially, if F denotes the vector field induced by (S, i?)^’ B \ it holds that Fx(V) = 
Fx(V) for all X £ p(S) and V £ R> 0 that are constant on H. 
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Partition Refinement Algorithms for CRN Bisimulations 


We study a polynomial-time algorithm for the computation of the coarsest bisimulations that 
refine an arbitrary input partition. We start introducing two auxiliary equivalence relations. 

► Definition 24 (Splitter equivalences). Let ( S,R ) be a CRN and H a partition over S. Then, 
we write A' Y if ([lj is fulfilled by (X, Y). Similarly, write X Y if (X, Y) satisfies (j2j. 

Algorithm [l] iteratively computes the coarsest forward or backward bisimulation (when 
X = F or x = B, respectively) that refines a given input partition of species of a CRN. 
Note that, contrary to CRN reduction algorithms, one (parametric) algorithm suffices for 
both bisimulations. Using the above splitter equivalences, at each iteration the blocks of the 
current partition S/~u are split in sub-blocks of '^-equivalent species 5/(~^ D ~n)- The 
algorithm terminates when no refinement is performed. 
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Algorithm 1 Template partition refinement algorithm for the construction of the coarsest 
CRN bisimulations that refine some given initial partition Q. 

Require: A CRN (S', R), a partition Q of S and % £ {F, B}. 

H <— Q 

while true do 

if %' = % then 
return H 
else 

Hi — W 

end if 
end while 


The freedom in choosing the initial partition Q is useful in both bisimulations. For FB it 
allows to single out species that are the “observables” of the CRN. These are the species for 
which the modeler is interested in obtaining distinct ODE solutions, information which would 
otherwise be lost if such species are found in larger equivalence classes. BB is lossless, hence 
this issue does not arise. Flowever BB requires the same initial conditions for equivalent 
species. In this case, an appropriate input partition may tell apart species for which it is 
known that the initial conditions are different. 

► Theorem 25 (Correctness). Given a CRN ( S,R) and a partition Q of S, Algorithm [I] 
calculates the coarsest forward and backward bisimulation that refines Q. In both cases, the 
number of steps needed is polynomial in the number of species and reactions. 

Note that, due to space constraints, we only focussed on the existence of a polynomial-time 
algorithm, and in the next section we provide numerical evidence of its scalability. The 
proof of this theorem gives a bound of 0(\R\ 2 ■ l^l 5 ) on the number of steps. Tighter bounds 
could be obtained by extending classical partition refinement approaches available for labeled 
transitions systems mm to CRNs, which is however the subject of future work. 
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Evaluation 


We now evaluate FB and BB. We first study their effectiveness in reducing the ODEs of a 
number of biochemical models from the literature given in the .net format of BioNetGen [5], 
version 2.2.5-stable. Using selected models we discuss how FB and BB relate with each 
other, and provide a biological interpretation of the aggregations. Finally, we compare them 
against /c’s fragmentation. All experiments are replicable using a prototype available at 
http://sysma.imtlucca.it/crnreducer/ 


Numerical results. Table [T] lists our case studies: four synthetic benchmarks to obtain 
combinatorially larger CRNs by varying the number of phosphorylation sites (M1-M4) |521 : 
a model of pheromone signaling (M5, |.‘dlj ): two signaling pathways through the Fee complex 
(M6-M7, HUES) ; two models of enzyme activation (M8-M9, [2]); a model of a tumor 
suppressor protein (M10, [3]); a model of tyrosine phosphorylation and adaptor protein 
binding (Mil, [T31 III] ): a MAPK model (M12, [21]); and an influence network (M13, [IT]). 

Headers \R\ and |5| give the number of reactions and species of the CRN (and of its 
reductions), respectively. The reduction times (Red.) account also for the computation of the 
quotient CRNs. The speed-up is the ratio between the time to solve the ODEs of the original 
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Original model 

Forward reduction 


Backward reduction 

Id 

\R\ 

|S| 

Red.(s) 

\R\ 

|S1 

Speed-up 

Red.(s) 

\R\ 

\s\ 

Speed-up 

Ml 

3538944 

262146 

4.61E+4 

990 

222 

— 

7.65E+4 

2708 

222 

— 

M2 

786432 

65538 

1.92E+3 

720 

167 

— 

3.68E+3 

1950 

167 

— 

M3 

172032 

16386 

8.15E+1 

504 

122 

1.16E+3 

1.77E+2 

1348 

122 

5.34E+2 

M4 

48 

18 

1.00E-3 

24 

12 

1.00E+0 

2.00E-3 

45 

12 

1.00E+0 

M5 

194054 

14531 

3.72E+1 

142165 

10855 

1.03E+0 

1.32E+3 

93033 

6634 

1.03E+0 

M6 

187468 

10734 

3.07E+1 

57508 

3744 

1.92E+1 

2.71E+2 

144473 

5575 

3.53E+0 

M7 

32776 

2506 

1.26E+0 

16481 

1281 

6.23E+0 

1.66E+1 

32776 

2506 

X 

M8 

41233 

2562 

1.12E+0 

33075 

1897 

1.12E+0 

1.89E+1 

41233 

2562 

X 

M9 

5033 

471 

1.91E-1 

4068 

345 

1.04E+0 

4.35E-1 

5033 

471 

X 

M10 

5797 

796 

1.61E-1 

4210 

503 

1.47E+0 

7.37E-1 

5797 

796 

X 

Mil 

5832 

730 

3.89E-1 

1296 

217 

1.32E+1 

6.00E-1 

2434 

217 

7.55E+0 

M12 

487 

85 

2.00E-3 

264 

56 

1.88E+0 

6.00E-3 

426 

56 

1.31E+0 

M13 

24 

18 

1.20E-2 

24 

18 

X 

7.00E-3 

6 

3 

1.00E+0 


Table 1 Forward and backward reductions and corresponding speed-ups in ODE analysis. Speed¬ 
up entries “—” indicate that the original model could not be solved; entries “x” indicate that the 
coarsest bisimulation did not reduce the original model. 


CRN and that of the reduced one including the time to reduce the CRN. Measurements 
were taken on a 2.6 GHz Intel Core i5 with 4 GB of RAM. The time interval of the ODE 
solution was taken from the original papers; for M1-M4, where this data was not available, 
time point 50.0 was used as an estimate of steady state. The initial conditions for the ODEs 
were also taken from the original papers. The initial partition for FB was chosen to be the 
trivial one containing the singleton block {5} (i.e., no species was singled out). Instead, the 
initial partition for BB was chosen consistently with the ODE initial conditions; that is, two 
species may be equivalent only if they have the same initial conditions in the original CRN. 
This ensured that the backward reduced CRN was a lossless aggregation of the original CRN. 

We make three main observations: (i) FB and BB can reduce a significant number of 
models. In the two largest models of our case studies, Ml and M2, the bisimulations were 
able to provide a compact aggregated ODE system which could be straightforwardly analyzed, 
while the solutions of the original models did not terminate due to out-of-memory errors, 
consistently with m- (ii) FB and BB are not comparable in general. For instance, both 
reduce M5 to 10855 and 6634 species, respectively, while M6 is reduced to 3744 species by FB, 
and to 5574 by BB. Also, FB was able to reduce M7-M10, while BB did not aggregate. The 
influence network M13 shows the opposite; in fact, none of the influence networks presented 
in [TT] can be reduced up to FB (here we showed Ml3, which is the largest one from [TT]'). 
(iii) Models M1-M4 and M12 show that the intersection between FB and BB is nonempty. 


Biological interpretation. Models Ml and M2 enjoy significant reductions and ODE analysis 
speed-ups. Here we use them to explain that FB and BB are effective at aggregating species 
representing symmetric sites in a complex. For this, let us consider M4, chosen for space 
reasons. A typical equivalence class is for instance {E(s\l).S(pl~P,p2~U\l), E(s!l).S , (pl~ 
E/!l,p2~P)}. According to the syntax of the BioNetGen language, the CRN species are 
formed from basic molecules S and E. Molecule S has two binding sites (pi, and p2) which 
can be either in phosphorylated state (P) or not (17); E has one stateless binding site (s) which 
can bind to those of S to form a complex. The two sites of S have equivalent capabilities in 
terms of binding with other species or changing state. For instance, the above equivalence 
class contains two species composed by S and E, with E bound to the unphosphorylated 
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site of S (here the exclamation mark links the binding sites used to form the species). 
Models Ml and M2 exhibit a fast growth of the number of species due to a larger number of 
symmetric sites, requiring distinct species to track exactly which site exhibits a particular 
phosphorylation state. This form of symmetry has also been studied in [3] where the authors 
propose an approach to detect it directly at the k level. However, an experimental comparison 
could not be performed because [3] is not yet implemented. Although both bisimulations 
give the same equivalence classes in these cases, the reduced CRNs have different reactions, 
since FB provides the dynamics of the sums of equivalent species, while BB considers the 
distinct dynamics of representative species. Instead, aggregation of identical binding sites 
is supported by BioNetGen. This can be seen in models M6 and M7, since they both have 
Lig(l : l ), a ligand with two copies of site l. Intuitively, the rule 

Rec(a ) + Lig(l , l) —> Rec(a\l).Lig(l\l, l ) (3) 

gives rise to only one chemical complex in the underlying CRN, Rec(a\V).Lig{l\l, l). This 
represents the (forward and backward) canonical representative of a ligand bound to a single 
receptor Rec(a). To see this, let us rename the two sites and expand the rule appropriately: 

Rec(a) + Lig(l i, I 2 ) —► Rec(a\l).Lig(li\l, I 2 ) 

Rec(a) + Lig(li, I 2 ) —► Rec(a\l).Lig(li, foil) (4) 

Then, this underlying CRN will distinguish the two sites. However, applying either of our 
CRN bisimulations leads to the CRN for Equation ([ 3 ]). 

We remark that the original CRN sizes of M 6 and M7 already account for the aggreg¬ 
ations obtained with BioNetGen. Nevertheless, our CRN bisimulations allow for further 
(significant) reductions. For instance, part of the reductions for M 6 are due to the presence 
of Rec{a , &, < 71 , < 72 ), a molecule with symmetric sites g\ and 52 , similarly to those of M4. 

Symmetric sites are not the only property captured by our bisimulations. For instance in 
both M8 and M9 one of the FB equivalence classes is given by: 

{ J(k\l).R(x\l, on, l), J(h\l).L(r\\2, r 2 )-R(x\l, on, Z!2), 

J(/c!l).L(ri, r2!2).i?(a;!l, on, Z!2), 

J(fc!l).T(ri!2, r 2 !3).R(x!l, ?'~on, l\3).R(x, z~on, d2)}. 

A biological interpretation is that a species containing the molecule J behaves in the same 
way as long as it is bound to a molecule R having binding site i in state “on”. This is 
independent of whether R is further complexed with other molecules via its binding site Z; 
For instance, the first species models that R is only bound to J, while in the second and 
third species it is also bound to L. Finally, in M5, one of the BB equivalence classes is 

{Dig2(p!l).Stel2(digl,dig2,dna!l,mapk), Fus3(p!l).Stel2(digl,dig2,dna!l,mapk), 
Msg5(p!l).Stel2(digl,dig2,dna!l,mapk), Sst2(p!l).Stel2(digl,dig2,dna!l,mapk), 
Stel2(p!l).Stel2(digl,dig2,dna!l,mapk), Ste2(p!l).Stel2(digl,dig2,dna!l,mapk)}. 

It captures that genes Dig2, Fus3, Msg5, Sst2, Stel2, and Ste2, bind to the protein Stel2 
with equal rates. This yields equivalent dynamics for these Stel2-gene complexes, and all 
those formed by them which are equal up to the gene bound to Stel2. 

Experimental comparison with K-based reduction techniques. We now experimentally 
compare our CRN bisimulations and fragmentation in the case of rule-based biochemical 


12 


Forward and Backward Bisimulations for Chemical Reaction Networks 


models for which the underlying CRN can be fully enumerated. All models in |Tablc T] belong 
to this class; however, none of them was originally available in k, the only language that 
supports fragmentation. Thus, we performed a manual translation of a selection of the case 
studies from the BioNetGen language into k. [^] We found: 

h Models that can be reduced by CRN bisimulations but not by fragmentation. The n 
encoding of M12 (a case where only cosmetic syntactical changes are required) returned 
85 fragments, equal to the size of the CRN, while both FB and BB reduced to 56 species. 
The encodings of M6 and M7 necessitated expansions analogous to Equation Q because 
k does not currently support distinct sites with the same name. This led to bigger initial 
CRNs, for which fragmentation returned 58040 fragments for M6 and 10930 for M7. 
b Models that can be reduced by fragmentation but not by our bisimulations. The n model of 
early events of the EGF pathway in [B] is reduced from 356 species to 38 fragments PE 
while no aggregation is obtained with either FB or BB. 
h Models that can be reduced by both our bisimulations and fragmentation. The n encodings 
of models M1-M4 present different reductions than using either bisimulation, specifically 
38, 34, 30 and 10 fragments (versus 222, 167, 122, and 12 FB and BB equivalence 
classes, respectively). It can be shown that, in the latter examples, the reductions are 
complementary, in the sense that no two bisimilar species are included in the same 
fragment. While our bisimulations captured symmetric sites, fragments explain that the 
sites of S are independent , i.e., the state of a site does not affect the dynamics of the 
other. For instance, one of the fragments for model M4 is 

{S(pl~P,p2~P), S{pl~P,p2~U), E(s\l).S(pl~P,p2~U\l), F(s\l).S(pl~P,p2~Pll)} 
which essentially collects all species where the pi site of molecule S is phosphorylated. 
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Conclusion 


Forward and backward bisimulations are equivalence relations over the species of a chemical 
reaction network inducing a partition of the underlying mass-action system of ordinary 
differential equations. An experimental evaluation has demonstrated their usefulness by 
showing their complementarity as well as significant model reductions in a number of 
biochemical models available in the literature. This has been supported by a prototype, 
which currently allows a ready-to-use tool-chain with BioNetGen, a state-of-the-art tool. 

Ongoing work is studying stochastic counterparts of both forward and backward bisimula¬ 
tions, to obtain model reductions when the semantics of chemical reaction networks based on 
continuous-time Markov chains is considered. Also, we plan to investigate the applicability 
of our bisimulations in other model repositories, e.g., those using the well-known SBML 
interchange format (http: //sbml. org). 
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APPENDIX 

A.l Forward CRN Bisimulation 

The following auxiliary lemma will be needed in the proofs of Propositions [To] and [TB] 

► Lemma 26. Let I be an index set and let 7 Zi denote equivalence relations on S such that 
S/lZi is a refinement of some partition Q of S. Then, TZ = (U is such that S/1Z is a 
refinement of Q. 

Proof. We first note that TZ is an equivalence relation over S, as it is the transitive closure of 
the union of equivalence relations over S. For i £ I, set Hi = S/IZi and H = S/1Z. Then, for 
any (j/o, y{) £ TZ, there exist xq, ... ,Xk £ S such that xoRi 0 xiRi 1 ... Ri k _ 1 Xk with yo = Xo, 
yi = Xk and ij £ I for all 0 < j < k — 1. Moreover, Xq £ G for some (unique) G £ Q. We 
show that Xq, ... ,Xk £ G by induction. Since the base case j = 0 is trivial, let us consider 
the induction step j — 1 —> j. Then, Xj-iTZi :j _ 1 Xj implies the existence of some H £ Hi j _ 1 
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such that Xj-i,Xj £ H. Let Gj-\ £ Q be such that H C Gj-\. Since Xj-i £ G by induction 
hypothesis and Xj-\ £ Gj-i, it holds that G ft Gj -1 ^ 0. Since Q is a partition, this implies 
that G = Gj- 1 , yielding in turn Xj £ G. ◄ 

Proof of Proposition |10[ We first note that TZ is an equivalence relation over S, as it is the 
transitive closure of the union of equivalence relations over S. For i £ I, set H, = S/IZi 
and H = S/1Z. Note that, for all i £ /, any H £ H is a union of blocks of Hi. For any 
(yoiVi) €E H, there exist xq ,...,Xk £ S such that XqR 1q x-[ R tl ■ ■ ■ Ri k _ 1 Xk with yo = Xq, 
yi = Xk and ij £ I for all 0 < j < k — 1. Noting that XjRj.Xj^ for all 0 < j < k — 1, 
we infer that err [xj,p\ = crr[xj + i,p] and pr[xj, p, H] = pr[xj +1 , p,H] for all p £ MS{S) 
and H £ H. Thus, we infer that crr[a; 0 ,p] = err [xi, p] = ... = cvv[xk-i 7 p\ = crr[x fe ,p] 
and pr [x 0 ,p,H] = pr[^i,p, H) = ... = pr[:Efc_i, p, H] = pr[a;fc, p , H], The remainder of the 
claim follows then from Lemma [26] ◄ 

We now provide a sort of commutative property for err and pr for the case of singleton 
p, used to prove Theorem 13 

► Proposition 27. Let ( S,R ) be a CRN and X, Y £ S. Then it holds 


crr(X, Y) = crr(F, A') 
pr(X,Y,Z) = pr(Y,X,Z) 


(5) 

( 6 ) 


Proof. The statement follows from Definition [7] since ([5]) and ([6]) can be rewritten, resp., as: 
flX,F|}(A). Y, «={|A, 11(F). Yi a 

X+Y-^tt£R X-\-Y — —^7r£R 

{|A,F|}(X). Y *-n{Z) = i\X,Y\ s {Y)- Y 

X+Y-^-w&R X+Y — —£ttGR 


Proof of Theorem [TT] 

For ( S , R) a CRN, H C S, and V a concentration function for (5, R), we use Vh for ^ Y eH Vx, 
cmdX H otY h to denote any element of H. We first separate the influence exerted by each 
reaction to the concentration of a species in a negative one (the depletion rate), and in a 
positive one (the accretion rate). 

► Remark. Let (S,R) be a CRN, V a concentration function for ( S,R ), and X £ S. We 
denote the depletion and accretion rates of X due to a reaction of R as, respectively, 

Depl{p TT, X, V) = p(X) • a • Vy (Y) 

Yes 

Accr(p -1 7r, X, V ) = x(X) ■ a • ^ 

Yes 


We can then reformulate the definition of vector field of (S,R) given in Section 2.1 
F x(V)= Y (Accr(p-^Tr,X,V)-Depl{p-^ir,X,V)) 


as: 


}tt£R 


We now state that the aggregated accretion and depletion rates of an equivalence class of 
an FB can be written in terms of the aggregate concentrations. 
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► Proposition 28 (Aggregated depletion and accretion rate). Let ( S,R ) be a CRN, V a 
concentration function for (S, R ), and R an FB. Then, for any H £ R 

E E Depl(p -^4 7r, V) X) and E E Accr(p -^4 7r, V, X) 

xeH p ^ neR xeH p ^ eR 

can be written in terms of the aggregated concentrations Vh 1 ■ ■ ■ Vh„ ■ 

Proof. We start addressing the depletion rate case. We have 

E E Depl(p^n,V,X)= E E P(X) ■ a ■ H V y 

XeH a Yaffle \ a<z(q\ <* yc.<? 


MY) 

'Y 


p - Yn£R 


xeH P eMS(S) p ^ neR 

e e n w 


yes 


p(Y) 


E PW 


XEH p eMS(S)YeS p -^ neR 

Given that we have unary or binary reactions only, we can rewrite ([T]) as: 

EE^- E m(x)-a+ 

X & HY£S y ^ neR 

E E Vv ' V y E $YXUX)-a 

XeHY+Y'EMS(S) Y+y ,^ veR 

Now, ([8]) can be easily rewritten in terms of the aggregated variables only: 

E E V v E ww •«= E v * E «= E ^ • crr t x ’ 0 i> 

xeHYes Y ^ neR xeh x ^ eR xeh 

from which by the condition on err of Definition § we obtain crr[A' ff , 0] • V R . 
Instead, as regards © we have: 

E E V r- V Y E {\Y,Y'\}(X).a = 

X£H Y+Y'eMS(S) y +y'-^-ker 

E E V x-V Y E = 

xehx+yems(S) x+Y ^ neR 

E y x E Vy ■ crr[X, Y] = (by the condition on err of Definition^ 

XEH X+YeMS(S) 

E y x E V Y ■ cvv[X H ,Y] = 

XEH X+YeMS(S) 

E Vx E E Vy • crr[X H ,y] = (by Proposition [gjp 

XEH HEHYEH 

’E Vx E E Vy ■ crr[y. X H ] = (by the condition on err of Definition^ 

XEH HEHYEH 


E Vx E crr[y*,X ff ] £ Vy = 

XEH h gR yeH 

E Vx E err \Y”,X h ]-V s = 

Yc TJ rr, 


XEH 

E Vh -cvt[Y S ,X h ] Y y x = 
heh xeh 

E V H ■ crr[yy X H ] ■V h = V h Y V h- cvv[X H , Y»] 
heh heh 


(7) 

( 8 ) 
(9) 
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closing the case. 

We now address the accretion rate case. We have 


£ £ Accr(p-^Tr,V,X)= 


p(y) 

Y 


xeH 


P ->7T(E R 


xeH peMS(S) 


p— yneR 
r p(Y) 


YeS 


= £ £ U v y £ *(*)■« ( 10 ) 


XeH p eMS(S)YeS 


p— yneR 


Given that we have unary or binary reactions only, we can rewrite (101 as: 


£ £ vy £ <*) 


XeHYeS 


Y - >-n -efl 

+ £ £ V Y ■ V Y > E »(*) ■« 

XCH y+y. e A,5(£) y+Y'^->x£R 

Now, ( |II] ) can be easily rewritten in terms of the aggregated variables only: 

£ IV £ £ n(X) ■ a = Y y Y ■ P r[Y,<h,H] = £ £ W- pr[Y,0,77] , 


( 11 ) 

( 12 ) 


res xeH 


Y - >k£R 


Yes 


HeUYeH 


from which, by the condition on pr of Definition E we obtain Y^PeU P r [l£ 77] • V R 
Instead, as regards we have: 

£ V Y -Vy'Y £ ’«= 

y+y'£MS{S) xen Y+Y ,^ eR 

£ v y £ £ < x ) • a + 

Y+YeMS{S) xeH Y+Y ^ GR 

£ Vy-Vy'Y £ 7f(X) ■ a = (see below) 


Y+Y'eMS(S) s.t. YjYf 


xeH 


Y+Y 1 - >k£R 


o£^£ 2 £ <*)■<*+ 


Yes xeH 


Y+Y - >ir£R 


o£^ £ v Y > y £ *(*)■« = 


(13) 

(14) 


YeS Y'es s.t. . YjtY 1 xeH 


Y+Y' - eR 


\Y v y- y ’ h ] + \J 2 y Y £ Vy ■ P r I y ’ y '> n] = 

z Yes z Yes Y'es s.t. y^y' 

2 £^v £ V Y '■ pr[Y,Y', H] = 
z Yes Y'es 

\ £ £ Vy £ Vy ■ pr[Y, Y', H) = (by cond on pr of Definition 


HeUYeH Y'es 


Y V u £ Vypr[Y»,Y',H} = 


Hen Y'es 


£ V R £ £ Vy' • pr [Y^,Y',H] = (by Proposition [D?]) 


Heu HeUY'eH 
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Vy ■ pr [Y 1 ,Y B , H\ = (by conct on pr of Definition HP 

HenY’eH 

\ Y y H E pr [Y S ,Y 6 ,H\ Y Vy’ = \ E V H E V s -pr[Y S ,Y B ,H] 

Hen H<z-H Y'eH Hen h<zu 


Where (13) is multiplied by 2/2, while p4| 
considered twice (e.g., Y + Y' and Y' + Y). 


is divided by two because each multiset is 

◄ 


Finally we can provide the proof of Theorem El 


Proof of Theorem 111! The theorem can be proved by showing that, for each H £ R, the sum 
of the ODEs of the species in H can be expressed in terms of the aggregated concentrations 
of the blocks of R only. Such sum is: 


E E Accr(p 7r, V, X) — E E Depl(p 7r, V, X) 


xeH 


xeH 


p YneR p—-^neR 

Thus, the claim follows from Proposition [28] 


◄ 


Proof of Theorem [20l 

Given a CRN (5, R) and an FB R, we hereby provide the technical results relating the 
%-lumped ODEs of (S, R) and the ODEs of its (R, E)-reduction. In the following we use 
X H to denote the canonical representative of the species in a block H £ R. In addition, 
for X S S and 7r g XiS(S), we may use px for p(X) and p n for p(tt). Finally, given a 
concentration function V for ( S,R ), we use yZH,F) f or the corresponding concentration 
function for (5, R)^ H ' F \ having a component = Vh per block H g %. 

We start providing a proposition used in the proof of Theorem [20] 

► Proposition 29. Let ( S,R ) be a CRN, 1Z an FB, R = S/1Z, and p its choice function. Let 
V be a concentration function for (S, R). Then, for any H £ R we have 


E E Accr(p -^4 tt,X,V) = Y Accr{p-^Ap(ir),X H ,V (H ’ F) ) 


xeH 


p —-^7r eR 


p T7TGH 

P=p(p) 


E E Depl(p^ir,X,V)= Y Depl(p ^ p{-k),X h M H ’ F) ) 


(15) 


(16) 


xeH 


p - YneR 


p - r-n-eR 

p=p(p) 


Proof. We only address the Acer case, as the Depl one is similar (actually, the latter is 
simpler). By the proof of Proposition 28 we obtain that we can rewrite the left-hand side of 
Equation © as 


_E P r [Y S ,<6,H}-V f j+ (17) 

Hen 

Hen ~ 


Hen 


( 18 ) 
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Instead, we can rewrite the right-hand side of Equation (151 as 

V ■ a ■ v^ r > - 


Y - fir£R 

Y=y,{Y) 


E M*") • a • V^ • v™ 


(19) 

( 20 ) 


y ± +y 2 -K ea 

i'i=A‘(> , i)Al r a=M(i r 2) 


We close the proof showing that Equation © is equal to Equation ( |T9) and that Equation ( |l8| ) 
is equal to Equation (20). 

We focus on the first equality, which easily follows noting that Equation (19) can be 
rewritten as 

E V H ,F) E m* H ) • a = E ^E E • a =Y, v fr #] 

Hen „ Hen Z ^ H ,, „ Hen 


Y - ?tt£R 

Y=X» 


•tt£R 

Y=XH 


We now consider the case “Equation (18) equals Equation (20)”. We can rewrite Equa¬ 
tion |20l) as follows: 


E M**) • «■ • yf F) + 


Y 1 +Y 1 - f"n£R 

Yi=#i(Vi) 


E t \?H \ xAH,F) rr(H,F) 

AM* )■ a -■ Vy 2 


( 21 ) 

( 22 ) 


n + ^2- >n€R,Yl?Y 2 

r 1 =A»(n)AF 2 =M(’bi) 


Now, Equation ( |21[ ) can be rewritten as 

E v h■ Vh e m^)-« = 


Hen 


Y 1 +Y 1 - 

Yi=*# 


_^ ^ ^ 2 

E ' V H E E 7r (^) ' a = (multiplying by -) 

Hen Z ^ H a , 

trl Vl+Vi- 

Y^X 6 

\ E ^ ^ E 2 E •«= 


Hen 


zeH 


Yi+Yi - fir eR 

Y 1= XH 


E V 6 -V B -pr[X 6 ,X A ,H\ 


(23) 


Hen 


Instead, Equation (221 can be rewritten as, where we divide Equation (241 by 2 because 
each multiset is considered twice (e.g., Y\ + Y 2 and Y 2 + Y\): 

\Y. V « E v s E = 

H&n Hen,H^H y 1 + v 2 -^ 6 i^ 

Y 1 =X B AY 2 =X H 

!E 7 s E y gE E ^• a = 


HeU,Hy£H 


zeH 


Y X +Y 2 -?-7refi 


(24) 











20 


Forward and Backward Bisimulations for Chemical Reaction Networks 


X E V ff E v s -pr[x*,x s ,ir\ 


(25) 


H&u 


By summing Equation (23) and Equation (251, we finally rewrite Equation (201 as 


IT i V h- v h-P^[X 6 ,X s ,H] + IY i V h E V s -pr[X S ,X S ,H} = 


Yen 




E V fr E Vg-pr [X”,X S ,H\ 


H&H hg-h 


This closes the case, and thus the proof is complete. 

We now provide the proof of Theorem [20] 

Theorem 1201 We first address the correctness of the reduction. For any H £ W we have 

E p xiv) = 53 E Accr (p v ) - E E De p l (p v ) 


xeH 


xeH 


>7 reR 


xeH 


YneR 


◄ 


As regards F, by Definition [l8j for any H GR we obtain 
F^(1/ (w ’ f) )= ^ Accr(p^p( n),X H M n - F) )- 


p 7 tvER 

p=m(p) 


E Depl(p^p(n),X H ,V^ F) ) 


p 7 tvER 

p=p(p) 


We close the proof by showing that 


E E Accr(p-?A ir,X,V)= Y Accr(pp{n), X H ,V {n ' F) ) 


xeH 


p —^7r eR 


p 7 tvER 

p=p(p) 


Y E D epl(p- 4tt,X,V)= 53 Depl(p^p(n),X H X H ’ F) ) 



(26) 

p(n),X H ,V^ F) ) 

(27) 


xeH 


p—Yn eR 


p 7 tvER 

p=p(p) 


Both Equations (26 ) and (271 follow from Proposition 29 


We now address the complexity of the reduction, showing that (S, R)^' H ’ F " > can be 
performed in 0(|F| • \S\ • (log(|5|) + log(|i?.|))) time. 

Steps (FI) and (F2) of Definition [l8| require to iterate (once) the reactions, (0(|i?|)). In 
particular, for each reaction we have to in turn iterate its reagents and products to perform 
(FI) and (F2), respectively. This requires 0(151) time. Finally, in order to efficiently perform 
(F3) we assume that the reagents and products are stored as an ordered list, and thus we 
have to sort the obtained canonized products (0(|5| • log(|5|))). To sum up, steps (FI) and 
(F2) require 0(|i?| ■ |S| ■ log(|5|)) time. 

Instead, step (F3) can be computed by first sorting the reactions obtained from (FI) and 
(F2), requiring 0(\R\ ■ log(|i?|) ■ |S|), where the |S| factor comes from the fact that in order to 
compare two reactions it is necessary to scan (once) their reagents and products. Then, (F3) 
is completed by iterating (once) the reactions to actually collapse them (<D(|i?| ■ |SQ). ◄ 
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A.2 Backward CRN Bisimulation 

Proof of Theorem [6j We first prove the if direction. Let p : S —> S be some choice function 
of F, set S := p(S') and define G^{V) := F%(V o p) for any V £ and X £ S. Further, 
let V denote the unique ODE solution of = G(V(t )) subject to some given initial 

condition F(0). Then, for all X £ S , it holds that 


°»)x = (^W) MY) = = F KX) { E(i) op) = F x {V{t) op) 

Thus, 1i-f V(t) o )jl is the unique solution of the ODE system ^V(t) = F(V(t)) subject to 
F(0) o fi. Since t i-> V(t) o p is constant on F, the proof of the if direction is complete. We 
now turn to the proof of the only-if direction. For this, let us assume that F is such that, for 
any F(0) £ R> 0 that is constant on F, the underlying solution of V = F(V) is constant on 
F as well. Fix arbitrary H £ F, X,Y £ H and V(0) £ R> 0 with V(0) constant on F. Since 
the solution of V = F(V) is differentiable, there exist S > 0 and functions r x ,r Y such that 
lim^o i r x{h) = lim^o \ r v{h) = 0 and 

V x (h) = Ky(0) + F y (K(0)) • h + r x (h) V Y {h) = W(0) + F Y (V( 0)) • h + r Y (h) 


for all 0 < h < S. Since V(0) is constant on F and X,Y £ if for some H £ F, this yields 
lim^o i(V x {h) - V Y (h)) = F X {V{ 0)) - FV(V(0)). Noting that V x (h) = V Y (h) for all 
0 < h < S because V(0) is constant on F, we thus infer Fx(F(0)) — Fy(V(0)) = 0. ◄ 

Proof of Proposition |16[ We first note that 7 Z is an equivalence relation over S, as it is the 
transitive closure of the union of equivalence relations over S. For i £ /, set F t = S/IZi 
and F = S/1Z. Note that, for all i £ I, any FI £ F is a union of blocks of Fi. From 
this, in turn, it is easy to see that any M £ {p \ p it £ R}/~ n can be written 
as a union of blocks of {p \ p tt £ R}/k,- h .. Observe that, for any (yo,yi) £ 7£, 
there exist Xo,-..,Xk £ S such that XoR^XiR^ ... Ri k l xu with yo = xq, y i = Xk and 
ij £ I for all 0 < j < k — 1. Noting that XjRi Xj+\ for all 0 < j < k — 1, we infer 
that fr [xj,M\ = fr^+i, M] for all M £ {p \ p it £ . Thus, we infer that 

fr [xq,M\ = ir[x\,M] = ... = fr[xfe_i, M] = fr [xk 1 M\ for all M £ {p \ p -^4 it £ 

The remainder of the claim follows then from Lemma |26j ◄ 

Proof of Theorem [T7J Define \p\v ■= Jives V X ' X> and set Q := {p \ p -^4 it £ 

Fix some arbitrary Xi,Xj £ FI and FI £ F and note that 


Fx k (V) 


55 a(it(X k ) - p(X k ))lpj v = 


55 55 ir ( X k,P)Mv 

[po]eQ pe[p 0 ] 


( 55 fr ( X fc>P))lPoJv 

[po]eQ pe[po] 

c(JC fc ,[po]) 


whenever V is constant on F. Observe also that the function V *-> F Xk {V), where V is 
constant on F, defines a polynomial in \Q\ variables with the monomials { c(Xk , [po]) ■ [po]v I 
[po] £ Q) ■ At last, recall that the multi dimensional version of Taylor’s theorem implies 
that two real polynomials are equivalent if and only if they have the same monomials. Thus, 
F Xj (V) = F x (V) for all V that are constant on F if and only if c(W, [p 0 ]) = c(Xj, [p 0 ]) for 
all [po] £ Q. ◄ 
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Proof of Theorem 1231 By encoding p : S —> S, p and 7r as arrays of length |S|, it is easy to 
see that the first operation needs at most 0(\R\ ■ |5|) steps. For the second operation, note 
that the renaming of species according to f .z can be done in again in 0(|R| ■ |S'!). However, since 
elements of MS(S) are stored as ordered lists to allow for performant processing, the second 
operation needs C>(|R| ■ |<Sj • log(|S'|)). To accomplish the third operation, instead, we first 
sort the reactions with respect to the lexicographical order which takes 0(\R\ • log(|R|) ■ IS 1 !). 
Afterwards, the rates of the reactions that coincide in reactants and products can be summed 
in 0(\R\ • |S|). 

We now turn to the correctness of algorithm. Let ( X , V ) i—>■ G l x (y) denote the vector field 
that arises from R after applying (Bl), ..., (Bz), with G° being F itself. We next prove that 
G' Xk (Vons) = Fx k {Vofj, s ) for all V £ R> 0 , Xj~ £ S and i £ {1,2}. For this, let us first apply 
the reaction changes p —> a i t ^ p —> a tt of (Bl). Then, G Xk (V o p) = G\ (V o p) because of 
n(Xk) — p(Xk) = 7r(Xfc) — p{Xk). Let us now consider a reaction changes p —> a n i->- p —> a tt 
of (B2). Then, G\JV o p) = G\ k (V o p) because ° A*) x X) = E[xes(^ ° 

and 7r(Xfc) — p{Xk) = n(Xk) — p(Xj~). Since G 2 Xk {Y ° p) = G 3 Xk (V o p) is trivially true, we 
infer the claim. ◄ 

A.3 Partition Refinement 

The following auxiliary results will be needed to prove the correctness of Algorithm [T] 

► Lemma 30. Given a CRN {S,R), let be two partitions of S such that Ri is a 

refinement ofH 2 . Then, the following holds true. 

1. Xi Xj implies Xi Xj. 

2. Xi ~}} i Xj implies Xi Xj. 

Proof. Let us assume that Xj Xj, i.e. it holds that err [Xi,p\ = err [Xj,p\ and 
pr[Xj, p, Hi] = pr[Xj, p, Hi\ for all p £ MS(S) and Hi £ Hi. Since any H 2 £ H 2 
can be written as a disjoint union of blocks from Hi, we thus infer that Xf Xj. 

Let us now assume that Xi Xj, i.e. it holds that tr[Xi,Mi] = ir[Xj,M\\ for all 
M\ £ {p | p tt £ R}/k U i . Since any M 2 £ {p \ p tt £ i?}/« W2 can be written as a 
disjoint union of blocks from {p \ p -M tt £ we deduce that Xj 

► Lemma 31. Let ( S,R ) be a CRN and H a partition of S. Then, the following holds. 

1. H is an FB if and only if Xi Xj for all Xi,Xj £ H and H £H. Moreover, it holds 

that H is an FB if and only ifH= fl ~«)- 

2. H is a BB if and only if Xj Xj for all Xi,Xj £ H and H £H. Moreover, it holds 

that H is a BB if and only ifH = fl ~^). 

Proof. The first parts of 1. and 2. are straightforward. The second parts, instead, follow 
from the corresponding first parts by observing that H = 5/(~^ D if and only if H is a 
refinement of S/ ◄ 

The following auxiliary results will be needed to establish polynomial complexity of 
Algorithm [TJ We implement p £ MS(S) as maps with keys in S. Moreover, we store a 
partition of a set A by means of a map that associates to each element of A a pointer to its 
representative. That is, each element of a partition block has a pointer to the representative 
of the block. 

► Lemma 32. Fix a CRN ( S,R ), pick A £ {S, R} and assume that deciding ai ~ a 2 for 
some equivalence relation ~ on A can be done in 0(|R| ei • |Sj 62 ) steps. Then, A/ ~ can be 
calculated in 0(|A| 2 • |R| ei • |R| e2 ) steps. 
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Algorithm 2 Algorithm to calculate the quotient A/ We assume that A is implemented 
as an array of objects {a[l],..., a[n]} where each object contains, among the actual data, a 
pointer that is initialized with zero at the beginning. 

Require: A set A and an equivalence relation ~ on A. 
for i = 1 to n do 

if a[i\.p != null then 
continue 
end if 

a[i].p 4— &a[i] 

for j = i + 1 to n do 

if a[j].p == null && a[j\ ~ a[i\ then 
a[j].p 4- ka[i\ 

end if 
end for 
end for 


Proof. It can be easily seen that Algorithm [ 2 ] calculates A/ ~ in 0(|A| 2 • |R| ei • |^1 62 ). ◄ 

► Lemma 33. For a CRN ( S,R) and a partition R of S, deciding X Y can be done in 

o(\r\ 2 -\s\ 2 ). 

Proof. We first note that, for a given p £ MS{S ), deciding err [A, p] = crr[Y, p\ can be 
done in 0(\R\ • |5|) because the comparison of two pi,p 2 £ AAS(S) needs 0(|S|). Similarly, 
for given p £ A4S(S) and Z £ S, the calculation of pr[X, p, Z) can be done in 0(\R\ • IS 1 !). 
Note also that T>(Z) := {p | 3 Z + p -^4 n £ R} can be calculated in 0(\R\ ■ |Sj) steps. We 
are now in a position to infer the claim. For this, note that err [X,p\ = err[Y, p] for all 
p £ MS(S) if crr[X, p] = crr[Y, p] for all p £ V(X)UD(Y). Consequently, deciding whether 
crr[X, p] = crr[Y, p] for all p £ MS(S) can be done in 0(\R\ 2 ■ |Sj) steps. Similarly, we note 
that pr[X, p , H] = pr[A, p, H] for all p £ MS(S) and H £ R if pr[X, p , H] = pr[Y, p, H] for 
all p £ V(X) U V{Y) and H £ "H. Thus, deciding whether pr[X, p, H] = pr[F, p, H] for all 
p £ MS(S) and H £R holds true can be done in 0(\R\ 2 ■ |Sj 2 ). ◄ 

► Lemma 34. For a CRN ( S,R ) and a partition R of S, deciding X Y can be done in 
0(|f?| 2 -|S|). 

Proof. Note that we can decide whether p a in 0(|»Sj). Thus, since 1Z = {p \ p 
7T £ R} can be calculated in 0(\R\ ■ IS 1 !) steps, Lemma [32] implies that Q = R,/ can be 
calculated in 0(\R\ 2 ■ |*S'|). Moreover, we note that, for a given p £ AAS(S), deciding whether 
fr[A, p] = fr [Y,p] can be done in 0(\R\ ■ |5|). Consequently, deciding whether X Y can 
be done in 0(\R\ 2 ■ |5|). ◄ 

We now are in a position to prove Theorem |25| 

Proof of Theorem 1251 For the proof of correctness, let us assume that R denotes the coarsest 
bisinrulation that refines Rq := Q and define Rk +1 := D ~n k ) f° r k > 0. Then, 

the sequence Ro,Ri,R 2 , • • • is such that 

1. R is a refinement of Rk 

2. Rk is a refinement of Rk- 1 

for all k > 1. We prove this by induction on k. 
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Listing 1 Encoding of M4 from BNGL into k , 

### SITE 1 ### 

#E (s) + S(pl~U) <-> E (s ! 1) . S (pi ~U! 1) 

E(s) , S(pl~U) <-> E(s!1),S(pl~U!1) 

#E (s ! 1) . S (p 1 ~U! 1) -> E(s) + S(pl~P) 

E(s!1),S(pl~U!1) -> E(s) , S(pl~P) 

#F(s) + S(pl~P) <-> F(s ! 1) . S (pl~P ! 1) 

F(s) , S(pl-P) <-> F(s ! 1) ,S(pl~P ! 1) 

#F (s ! 1) . S (pi ~P ! 1) -> F (s) + S(pl~U) 

F (s ! 1) , S(pi~P!1) -> FCs) , S(pl~U) 


### SITE 2 ### 

#E(s) + S(p2~U) <-> E(s ! 1) . S (p2~U! 1) 
E(s) , S(p2~U) <-> E(s ! 1) ,S(p2~U!1) 

#E (s ! 1) . S (p2~U! 1) -> E(s) + S(p2~P) 

E(s!1),S(p2~U!1) -> E(s) , S(p2~P) 

#F(s) + S(p2~P) <-> F(s ! 1) . S (p2~P ! 1) 
FCs) , S(p2 ~P) <-> FCs ! 1) ,S(p2~P!1) 

#F(s ! 1) . S(p2~P ! 1) -> FCs) + S(p2~U) 
FCs!1) ,SCp2~P!1) -> FCs) , S Cp2 ~U) 


k = 1: Since H is a refinement of Ho, Lemma 30 ensures the first claim. The second 
claim is trivial. 

k —> k + 1 : Thanks to the fact that % is a refinement of Hk by induction, Lemma 30 
ensures the first claim. The second claim is trivial. 

Since H is a refinement of any H/c, it holds that H = Hk whenever Hk is a bisimulation. 
Thanks to the fact that Hk is a refinement of Hk-i for all k > 1 and S is finite, we can fix the 
smallest k > 1 such that Hk = Hk- 1 - Since this implies Hk -i = Hk = S/{~u k -i ^ 

Lemma [31] yields the claim. We now turn to the complexity analysis. Note that since deciding 
X Y needs a constant number of steps, Lemma 33 and Lemma 34 imply that X 
can be decided in(9(|H| 2 -|S'| 2 ) steps. Thus, Lemma 32 ensures that S/ 
in 0{\R\ 2 • |S'! 4 ) steps. Noting that the algorithm makes at most |S| iterations, we conclude 
that the algorithm needs at most G(\R\ 2 • |5| 5 ) steps. ◄ 


x y 
n 1 

U can be calculated 


A.4 ^-encodings discussed in Section [6] 

M1-M4 


We start providing the K-encoding of model M4 of Table 1 whose original (BNGL) version 
has been taken from [32]. M4 is a special case in which almost no changes are necessary 
to convert it in k. The encoding is given in [Listing 1[ where we omit unnecessary details. 
Each k rule is preceded by the corresponding original BNGL rule (starting with #). The 
encodings for all other models of the same benchmark, M1-M3, are similar and thus we omit 
them here. They are available for download at http://sysma.imtlucca.it/crnreducer/ 


M6-M7 

As discussed in Section [6] models M6 and M7 are not directly encodable in k, as they contain 
the molecule Lig(l , l) having two identical binding sites l, which is forbidden in k. Hence, we 
encoded in k an expansion of the models similar to that of Equation Q. For presentation 
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reasons, we do not provide the /c-encodings of the expansions of M6 and M7, which however 
are available for download at http://sysma.imtlucca.it/crnreducer/ 


M12 


We conclude this appendix providing in Listing 2 the K-encoding of M12, whose BNGL 
version has been taken from m- As for M4, each k rule is preceded by the corresponding 


original BNGL rule. 
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Forward and Backward Bisimulations for Chemical Reaction Networks 


Listing 2 Encoding of M12 from BNGL to k 

#MAP3K (s,S~I) -> MAP3K(s,S~A) 

MAP3K(s,S~I) -> MAP3K(s,S~A) 

#MAP3K(s,S~A)+Scaff(map3k) <-> MAP3K(s!1,S~A) . Scaff(map3k!1) 

MAP3K(s,S~A),Scaff(map3k) <-> MAP3K(s!l,S~A),Scaff(map3k!1) 

#MAP3K(s!1,S~I).Scaff(map3k ! 1) -> MAP3K(s,S~I)+Scaff(map3k) 

MAP3K(s!1,S~I) ,Scaff(map3k ! 1) -> MAP3K(s,S~I) ,Scaff(map3k) 

#MAP2K (s ,R1~Y, R2~Y)+Scaff (map2k) <-> MAP2K(s ! 1, R1 ~Y, R2~Y) . Scaff (map2k ! 1) 

MAP2K(s,R1~Y,R2~Y),Scaff(map2k) <-> MAP2K(s!1,R1~Y,R2~Y),Scaff(map2k!1) 

#MAP2K (s , Rl~Yp , R2 ~ Y) +Scaff (map2k) <-> MAP2K(s ! 1, R1 ~Yp , R2~Y) . Scaff (map2k ! 1) 

MAP2K(s,Rl~Yp,R2~Y),Scaff(map2k) <-> MAP2K(s!l,Rl~Yp,R2~Y),Scaff(map2k!1) 

#MAP2K (s ,R1~Y, R2 ~Yp ) + Scaff (map2k) <-> MAP2K(s ! 1, R1 ~Y, R2~Yp ) . Scaff (map2k ! 1) 

MAP2K(s,Rl~Y,R2~Yp),Scaff(map2k) <-> MAP2K(s!1,R1~Y,R2~Yp),Scaff(map2k!1) 

#MAP2K(s ,Rl~Yp , R2 ~Yp ) + Scaff (map2k) <-> MAP2K(s !1, Rl~Yp , R2~Yp) . Scaff (map2k ! 1) 

MAP2K(s ,Rl~Yp ,R2~Yp) , Scaf f (map2k) <-> MAP2K(s!l,Rl~Yp,R2~Yp) , Scaf f (map2k ! 1) 

#MAPK(s ,R1~Y, R2~Y)+Scaff (mapk) <-> MAPK(s ! 1 , R1~Y, R2~Y) . Scaff (mapk ! 1) 

MAPK(s,R1~Y,R2~Y),Scaff(mapk) <-> MAPK(s!1,R1~Y,R2-Y),Scaff(mapk!1) 

#MAPK(s ,Rl~Yp , R2~Y)+Scaff (mapk) <-> MAPK(s ! 1, R1 ~Yp , R2~Y) . Scaff (mapk ! 1) 
MAPK(s,Rl~Yp,R2~Y),Scaff(mapk) <-> MAPK(s!1,R1~Yp,R2~Y),Scaff(mapk!1) 

#MAPK (s , R1 ~Y, R2~Yp)+Scaff (mapk) <-> MAPK(s ! 1, R1 ~Y, R2~ Yp ) . Scaff (mapk ! 1) 
MAPK(s,Rl~Y,R2~Yp),Scaff(mapk) <-> MAPK(s!1,Rl~Y,R2~Yp),Scaff(mapk!1) 

#MAPK(s !1 ,Rl~Yp ,R2~Yp) . Scaf f (mapk ! 1) -> MAPK (s , R1 ~Yp , R2~Yp ) + Scaf f (mapk) 

MAPK(s!1,Rl~Yp,R2~Yp),Scaff(mapk!1) -> MAPK(s,R1~Yp,R2~Yp) , Scaff(mapk) 

#MAP3K(s !1, S~A) . Scaff (map3k !1, map2k !2) . MAP2K(s !2,R1~Y) -> MAP3K(s ! 1, S~A) . Scaff (map3k ! 1, 
map2k !2) . MAP2K(s !2,Rl~Yp) 

MAP3K(s!1,S~A),Scaff(map3k!1,map2k!2),MAP2K(s!2,R1~Y) -> MAP3K(s!1,S~A),Scaff(map3k!1, 
map2k!2),MAP2K(s!2,Rl~Yp) 

#MAP3K(s !1, S~A) . Scaff (map3k !1, map2k !2) . MAP2K(s !2,R2~Y) -> MAP3K(s ! 1, S~A) . Scaff (map3k ! 1, 
map2k !2) . MAP2K(s !2,R2~Yp) 

MAP3K(s!1,S~A) ,Scaff(map3k!1,map2k!2) ,MAP2K(s!2,R2~Y) -> MAP3K(s!1,S~A) ,Scaff(map3k!1, 
map2k!2),MAP2K(s!2,R2~Yp) 

#MAPK(s ! 1, R1 ~Y) . Scaff (mapk ! 1, map2k ! 2) . MAP2K(s 12, Rl~Yp , R2~Yp) -> MAPK(s !l,Rl~Yp) . Scaff( 
mapk ! 1, map2k !2) . MAP2K(s !2,Rl~Yp , R2~Yp) 

MAPK(s!1,R1~Y) ,Scaff(mapk!1,map2k ! 2) ,MAP2K(s!2,R1~Yp,R2~Yp) -> MAPK(s!1,R1~Yp) ,Scaff( 
mapk ! 1 , map2k ! 2) , MAP2K ( s ! 2 , R1 - Yp , R2 ~ Yp) 

#MAPK(s ! 1 , R2~Y) . Scaff (mapk ! 1 , map2k ! 2) . MAP2K(s 12, Rl~Yp , R2~Yp) -> MAPK(s !l,R2~Yp) . Scaff( 
mapk ! 1, map2k !2) . MAP2K(s !2,Rl~Yp , R2~Yp) 

MAPK(s!1,R2~Y) ,Scaff(mapk!1,map2k!2) ,MAP2K(s!2,R1~Yp,R2~Yp) -> MAPK(s!1,R2~Yp) ,Scaff ( 
mapk!1,map2k!2),MAP2K(s!2,R1-Yp,R2~Yp) 

#MAP3K(S~A) -> MAP3K(S~I) 

MAP3K(S-A) -> MAP3K (S~I ) 

#MAP2K(Rl~Yp) -> MAP2K (R1 ~ Y) 

MAP2K(R1~Yp) -> MAP2K(R1~Y) 

#MAP2K (R2~Yp) -> MAP2K (R2~Y) 

MAP2K(R2 ~Yp) -> MAP2K(R2~Y) 

#MAPK (R1 ~Yp) -> MAPK(Rl-Y) 

MAPK(R1~Yp) -> MAPK(Rl-Y) 

#MAPK (R2~ Yp ) -> MAPK (R2 ~ Y) 

MAPK(R2~Yp) -> MAPK(R2~Y) 



